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We review progress towards direct simulation of quantum dynamics in many-body 
systems, using recently developed stochastic gauge techniques. We consider master 
equations, canonical ensemble calculations and reversible quantum dynamics are 
compared, as well the general question of strategies for choosing the gauge. 



1. INTRODUCTION 

In this paper, we will review progress toward the direct simulation of quan- 
tum dynamics in many-body systems. This is a central problem in modern 
theoretical physics, sometimes claimed to be inherently insoluble 1 on dig- 
ital computers, due to complexity issues. In particular, we will evaluate 
progress and results using stochastic gauge methods 2 on currently avail- 
able digital computers. 

This is a timely issue to the spectroscopist at the start of the twenty-first 
century, since spectroscopy now involves the study of increasingly complex 
systems - not just molecules, but large interacting many-body systems like 
dilute bosonic and fermionic quantum gases. The use of subtle dynamical 
experiments allows the exploration of many-body dynamics well beyond 
regimes where perturbation theory is useful, and experiments in regimes 
where mean-field theory is not applicable are increasingly common. 

The specific issues treated are: 

• Stochastic gauge - this is a unified simulation method which pro- 
vides a route to understanding a wide range of stochastic methods. 

• Master equations - this is the most general real time problem, which 
can include both reversible and/or dissipative elements. 

• Canonical simulations - needed if a finite temperature grand canon- 
ical ensemble is to be calculated, in thermal equilibrium. 

• Quantum dynamics - we illustrate this by a calculation of correlated 
atom pairs in coherent molecular conversion to atoms. 



1 
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1.1. The problem: complexity 

Before giving an account of these new techniques, we must consider what 
the problem is. As an example, consider modern BEC experiments, which 
may involve 10 100000 states, or around 10 6 equivalent qubits in Hilbert 
space. Such many-particle dynamical experiments are commonplace, but 
they do place severe demands on theory. 

Apart from various approximate methods, direct computation in a 
number-state basis is useful only for small numbers of particles; in gen- 
eral, it needs exponentially too much memory! The last problem led to a 
famous and somewhat pessimistic conclusion by Feynman 1 , summarized in 
the following statement: "Can a quantum system be probabilistically simu- 
lated by a classical universal computer? ..the answer is certainly, No!" 

1.2. Quantum Computers: too small, too costly? 

A possible solution to the complexity problem was suggested by Feynman. 
This is the proposal to develop quantum computers in which the logical 
entities are qubits - that is, two-state quantum systems - and the logical 
operations are dissipation-free unitary transformations on large numbers of 
qubits. This solves the quantum complexity problem at the software level 
- since qubits are extremely efficient at storing large Hilbert spaces - but 
there is still a hardware issue to be solved. 

Great progress has been made at NIST 3 (also reported at this confer- 
ence), in constructing ion trap devices that are able to carry out binary 
qubit logic operations. Despite this, even if the NIST ion traps are able to 
be scaled to the projected size of 10 4 ions, they will still be relatively small 
due to error-correction overheads, and also rather slow and expensive. This 
is illustrated in the following table. 



COMPUTER 


bits 


Hz 


Ops/cycle 


$/CPU 


Ops/s/$ 


DIGITAL(2000) 


10 10 


10 9 


10 2 


10 3 


10 8 


QC (2000) 


10 


10 3 


1 


10 6 


10" 3 


DIGITAL (2020) 


10 12 


10 10 


10 3 


10 


10 12 


QC (2020) 


10 2 


10 4 


10 


10 5 


1 



While the digital computer projections are industry-based 4 , the quan- 
tum computer projections are optimistic in assuming that all of the known 
problems will be solved on this timescale. 
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2. Phase-space representations 

There is a possible alternative: improved digital algorithms that may allow 
industry standard digital computers to be used for quantum simulations. 
The proposal treated here is to map quantum dynamics into equations for 
stochastic motion on a phase-space, which are then randomly sampled to 
obtain statistical averages of quantum observables. The procedure is similar 
to path-integral Monte-Carlo techniques 5 used for canonical ensembles. 

Historically, this method originated in the classical phase-space meth- 
ods: the Wigner 6 , Q 7 , and P 8 representations all have a classical phase- 
space dimension of (say) d dimensions. However, the quantum dynamics 
of systems of interacting particles cannot be sampled stochastically. More 
recently, methods of larger dimension using 'quantum' phase-spaces were in- 
troduced: the first was the Positive-P representation 9 , using a generalized 
phase-space variable a with d complex dimensions. This allows dynamical 
problems to be mapped into stochastic equations, given some restrictions 
on the tails of the distribution. 

The positive-P technique has been used widely in quantum optics, and 
led to the prediction of quantum soliton squeezing 10 , which was verified 
experimentally 11 . While the first calculations used linearization, these were 
quickly extended 10 to the first real-time quantum field simulations with 10 9 
particles in 10 3 modes, corresponding to a Hilbert space dimension of more 
than 10 4 qubits. 

More recently, these methods were extended to many other nonlinear 
quantum systems, including full two and three dimensional calculations, 
typically with damping present. An example of this is the first a-priori 
quantum simulation 12 of the formation of a BEC via evaporative cooling. 
With 2 x 10 4 particles in 3 x 10 4 modes, this corresponds to over 10 5 qubits. 

Despite its success, the basic positive-P technique does have limitations. 
The stochastic equations have stability problems if there is not sufficient 
damping in the system, leading to uncontrollable growth in sampling error 
after a certain simulation time. Systematic 'boundary term' errors 13 can 
arise for certain nonlinear systems. 

Both of these problems can be solved by use of stochastic gauges 2 . 
This extension of the positive-P technique introduces an extra phase-space 
variable 0, which behaves as a weighting function, and allows the intro- 
duction of arbitrary 'drift gauges' g(a), which can be used to stabilize 
stochastic paths without affecting the physical ensemble averages. In addi- 
tion, there are 'diffusion gauges' g d (a), which exploit a freedom in choice 



February 2, 2008 4:6 WSPC/Trim Size: 9in x 6in for Proceedings 



GAUGEICOLSpapcr2 



4 

of noise terms to reduce the sampling error. The stochastic gauge method 
is a very general technique, which also unifies some earlier methods ' 15 . 
Its flexibility enables it to be tailored to different applications. 

2.1. Method in outline 

A general quantum calculation in real or imaginary time (for thermal equi- 
librium) can be written as a Liouville equation for the density operator: 

dp(t)/dt = L[p(t)}- (1) 

To solve this with stochastic gauges, we first expand the density operator 
over some suitable operator basis A( A ): 

p(t) = J P{X ,t)A(^)d 2{d+1) ^ , (2) 

which defines a (d + l)-dimcnsional complex phase-space: A = (0,a). 
For a given basis, there are identities which allow us to write the Liouville 
operator equation as 

dp(t)/dt = J P{X,t)C A H~\)d 2{d+1 '>^, (3) 

where La is a linear differential operator which can include arbitrary gauge 
functions. If there are no derivatives higher than second order, this equation 
can be transformed to a Fokker-Planck equation for P, provided the gauges 
are chosen to eliminate any boundary terms that may otherwise arise from 
the partial integration. This can always be transformed into the stochastic 
equations: 

dn/dt = Q [U + g • C] (4) 
da/dt = A + B(C- g), (5) 

where U and the vector A are determined by the form of the original 
Liouville equation. The drift gauges appear as the arbitrarily functions g, 
and diffusion gauges appear as the freedom that exists in choosing the noise 
matrix B. The noise terms £ are Gaussian white noises, with correlations: 

m)Q(t'))=8 t3 6(t-t'). (6) 

Equations (4-6) form the central result of the stochastic gauge method 
and can be used to solve a large class of quantum dynamical and thermal- 
equilibrium problems 2 . 
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In practice, the numerical implementation of these equations is simpli- 
fied by use of an extensible multi-dimensional simulator (XMDS) which 
generates efficient C++ computer code from a brief high-level problem de- 
scription. It also features automatic clustering with standard message- 
passing interface (MPI) routines, and is free under the open-source GNU 
public license (GPL) at the website www.xmds.org. 

There are three main types of problems which can be studied using 
these methods: 1) master equations of open quantum systems, which fea- 
ture damping as well as coherent evolution, 2) canonical ensembles, which 
involve the 'imaginary time' calculation of thermal equilibrium states and 
3) quantum dynamics of closed many-body systems. We now present an 
example of each type that has been solved by use of stochastic gauge meth- 
ods. 

3. Master Equation 

Accurate simulation of experimental systems usually requires taking into 
account damping and other coupling to the external environment, which 
can be achieved with a master equation. The positive-P method has been 
shown to suffer from boundary term systematic errors 13 in some nonlin- 
ear cases when there is small linear damping, for example in a coherently 
driven interferometer with two-boson damping. Introducing appropriate 
drift stochastic gauges removes these systematic errors. 

A typical master equation , when reduced to a single mode, can be 
written 

^ = [eat _ e * a , p\ + \{2a 2 P a) 2 - a) 2 a 2 p - pa) 2 a 2 ) , (7) 

with e the driving field amplitude. The performance of stochastic gauge 
simulations has been compared to positive-P and exact results in earlier 
work 2 . Due to space limitations, we refer the reader to this article, which 
shows that in all the cases studied, any systematic errors are removed and 
the useful simulation time is greatly extended. The stochastic gauge method 
then becomes the preferred approach when applied to a many-particle or 
many-mode system in which a direct number-state basis calculation is im- 
practical. 

4. Grand canonical ensembles 

We consider the equilibrium states of an interacting Bose gas, thermally and 
diffusively coupled to a reservoir. The calculation of spatial correlations in 
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a one (or more) -dimensional interacting Bose gas is a highly non-trivial 
task, and has been calculated for the first time using the stochastic gauge 
method 14 ' 16 . The positive-P method is inadequate to the task because of 
boundary term errors, and the inability to handle the Gibbs weight factors. 

The imaginary time evolution of grand canonical ensembles (with tem- 
perature T) can be modeled using an anti-commutator equation of form: 



dp u 
8t 



N 



dp, 

a7 



H,p u 



(8) 



Here, r = l/k B T, p is the chemical potential, N is particle number, H 
the Hamiltonian, and p u is un-normalized. One starts at very high tem- 
perature r = 0, where the ensemble is known, and evolves in r to get 
finite-temperature ensembles. 

We model the ID Bose gas with the Bose- Hubbard Hamiltonian (h = 1) 



H = 



EE^a^+E : 



(9) 



on a lattice labeled by i. The w, 3 are nonlocal couplings between modes, 
and ni are occupations at each point. 
One obtains the equations 
don 

— = -[m + igi - i(,i(T)\ai - 2_^("iy < i 



dr 



[rii + ig, - iQ t {T)]a t 



iOj/2, 



ft/2, 



E«[CiW+Ci(r)]. 



(10) 



Here the include w, 3 together with terms coming from the chemical 
potential, and rij = a,/3j, while Ci( T ) an( i d( T ) are an independent Gaussian 
noises of variance 5(t). We choose the stabilizing gauge gi = i(Rc[rij] — \m\) 
to remove boundary term errors, and reduce sampling error. 
The second-order correlation function 



gV(x) = (:n(0)n(x) :) / (n(0)) (n(x)) 



(11) 



is shown in Fig. 1 for a uniform one-dimensional Bose gas in the Tonks- 
Girardeau regime of incipient fermionization. 

It is seen that g^ (x) is well calculated, and shows that in this parameter 
regime, there is a preferred inter-atom distance of approximately half the 
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Figure 1. Momentum distribution n(k) and correlation function (x) of uniform 1-D 
Bose gas (solid lines). Here 7 = 10 = mg/p is the scaled interaction strength (with p 
the linear density), and T = lOT^, where = 2itp 2 /mkB is the quantum degeneracy 
temperature and £ is the healing length. ID is for an ideal gas (7 = and T = lOT^), 
while BL is for a Boltzmann gas (7 = 10, T/Tj — > 00). 

healing length. Anti-bunching {g^'(0) = 0.72 ±0.01) is also seen. There is 
complete agreement with all known finite temperature exact results 17 . 



5. Quantum Dynamics 

Here we consider the process of coherent dissociation of a BEC of molecular 
dimers into pairs of constituent atoms, via stimulated Raman transitions 
or Feshbach resonances 18 . The resulting effect is formation of quantum 
correlated twin atom-laser beams with squeezing in the particle number 
difference. This is the matter wave analog of optical parametric down- 
conversion producing the famous entangled photon pairs. 

The effective Hamiltonian interaction, taken for simplicity in one space 
dimension, is 



■ X(t) 



I 



dx 



t2 



(12) 



Here, ^>i(x,t) and ^> 2 ( x ,t) are the atomic and molecular field operators, 
respectively, x(t) is the atom- molecule coupling responsible for the conver- 
sion of molecules into atom pairs (and vice versa), and w is a detuning. The 
complete Hamiltonian also contains the usual kinetic energy terms, exter- 
nal trapping potentials, as well as usual quartic interaction terms describing 
atom-atom, atom-molecule, and molecule-molecule s-wave scattering. 

The quantum dynamical simulation of this process is carried out using 
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the +P representation. This is an example where interesting physical ef- 
fects can take place on short time scales so that the straightforward +P 
simulations give reliable results. 

For short dissociation time intervals such that the atomic self-interaction 
can be neglected due to low atomic density, and for large effective detun- 
ing such that the atom-molcculc s-wave interaction can be neglected too, 
the dynamics is simulated via the following set of stochastic equations (in 
appropriate dimcnsionlcss units): 

dij)i .d 2 i>\ 
Ih = 

~d7 



1 d 2 tfj 2 

2 oe 



(7 + id)tpi - 
- w(£,r)V>2 



(13) 



together with the corresponding "conjugate" equations for the i/^-fields. 
Here, £ and r are the dimcnsionlcss coordinate and time, k is the coupling 
proportional to \, t) accounts for the molecular trapping potential and 
u the molecule-molecule self-interaction, 5 is the dimensionless detuning 
parameter, 7 accounts for possible atomic losses, and 771,2 are the Gaussian 
delta-correlated noise terms. 
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Figure 2. Evolution of the atomic density (left panel) and the variance V versus time 
(right panel). The dissociation takes place in the time interval from r = to 8 X 10 — 4 , 
followed by free evolution of the atomic beams. 



A typical example demonstrating the formation of two counter- 
propagating atomic beams is given in Fig. 2. The correlations between 
the beams propagating to the left (-) and to the right (+) are quantified 



February 2, 2008 4:6 WSPC/Trim Size: 9in x 6in for Proceedings 



GAUGEICOLSpapcr2 



9 

via the normalized variance of fluctuations in the relative total number of 
particles, V = ([A(N- - N+)} 2 ^ / ((-#-) + (^+))- In a coherent (un- 
correlated) state, the variance is V = 1, while V < 1 implies squeezing 
of quantum fluctuations below the coherent level - which is due to strong 
quantum correlations between the atoms in the two beams. The example 
represented in Fig. 2 shows more than 93% squeezing, once the dissociation 
is switched off and the atomic beams separate spatially. 

The physical reason for the correlation and the squeezing is momentum 
conservation, which requires that each atom emitted with a (dimensionless) 
momentum q > be accompanied by a partner atom having q < 0. In order 
to conserve energy, this pairing only occurs for 5 < 0, which allows the 
potential energy in the molecule to be converted to atomic kinetic energy 
for selected modes with q values around qo — ±y/\S\. 

6. Strategies and future developments 

The stochastic gauge representation is a very general technique suited to 
a large class of problems. The secret of its success lies in the possibility 
of adapting the equations to the particular problem being studied. Al- 
ready much work has been done on optimizing the gauge for a variety of 
problems. Based on this, we conjecture that in the general, a successful 
drift gauge should a) guarantee inward asymptotic flow, b) generate an 
attractive manifold on which the gauge itself vanishes, and c) be used in 
conjunction with a diffusion gauge to minimize sampling error. Using these 
techniques we have been able to extend the useful time of purely unitary 
(lossless) quantum time-evolution by two orders of magnitude compared to 
the positive-P method, for anharmonic oscillator time-evolution with up to 
10 10 bosons. 

Another, complementary strategy is to optimize the basis set. A basis 
set that incorporates more of the physics of the problem being studied 
would lead to more efficient sampling. Thus, we have recently formulated 
a phase-space technique that uses a Gaussian basis 19 , which is also suited 
to generating a phase-space method for fermions. 

A third strategy is to optimize the stochastic sampling, employing such 
QCD-likc techniques as the Metropolis and genetic algorithms. So far, these 
techniques have only just started to be investigated. 

Perhaps the most important overall conclusion is that while ab-initio 
quantum dynamical simulations are difficult, they are not ruled out in 
principle. These techniques can serve as useful benchmarks in establish- 
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ing the validity of approximate methods. They are not limited to quan- 
tum problems, since they can also be used for kinetic, genetic or chemical 
master equations. But the most interesting development will be in new 
tests of quantum mechanics for systems that are both strongly entangled 
and macroscopic - and these are the types of quantum state whose time- 
evolution is not easily predicted using previous methods. 
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